Simulation study comparing rule-based and model-based adaptive designs
James Watson
January, 2020
Setup model parameters
source('Core_functions.R')
source('Simulation_functions.R')
source('plotting_functions.R')
source('wrapper_function_all_scenarios.R')
TEL = 0.95 # Target Efficacy Level
MTT = 0.05 # Maximum Tolerated Toxicity
N_max = 260 # Maximum number of patients recruited
max_increment = 1 # Maximum dose increment to doses previously unseen
N_trials = 2000
N_batch = 4
Randomisation_p_SOC = 0.2 # proportion randomised to the standard of care dose
starting_dose = 12 # starting dose in adaptive arm
SoC = 8 # Standard of Care
# function that solves for the beta value based on interpretable parameters
solve_beta = function(alpha_val, v_star, y_star){
beta_val = ( logit(y_star) - alpha_val ) / log2(v_star)
return(beta_val)
}
optimal_doses = array(NA,dim = 7)Priors
#******* Prior point estimates *********
Prior_MTD = 32; # prior estimate of the Maximum Tolerated Dose
Prior_alpha_tox = logit(1/1000) # prior estimate of the toxicity after 1 vial
Prior_beta_tox = solve_beta(alpha_val = Prior_alpha_tox, v_star = Prior_MTD, y_star = MTT)
#******* Prior uncertainty estimates *******
prior_model_params = list(beta_tox = Prior_beta_tox,
beta_tox_sd = .05,
alpha_tox=Prior_alpha_tox,
alpha_tox_sd = 2,
mu_antivenom = 8,
mu_antivenom_sd = 3,
sd_antivenom = 5,
sd_antivenom_sd = 2)Simulation 1
true_TED = 20 #
true_antivenom_mean = 10
true_antivenom_sd = (true_TED - true_antivenom_mean)/qnorm(.95)
true_alpha_tox = logit(1/500) # toxicity at 1 vial
true_MTD = 8 # simulation truth for the MTD
true_beta_tox = solve_beta(alpha_val = true_alpha_tox, v_star = true_MTD, y_star = MTT)
# well-specified simulation truth
model_params_true = list(beta_tox = true_beta_tox,
alpha_tox=true_alpha_tox,
mu_antivenom=true_antivenom_mean,
sd_antivenom=true_antivenom_sd)
optimal_doses[1] = min(true_MTD, true_TED)
tic()
Full_Simulation(model_params_true = model_params_true,
true_model = NULL,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 1',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'model_based',
starting_dose = starting_dose,
SoC = SoC,
use_SoC_data = T)## [1] "Simulation scenario 1_model_based_well_specified_All_data.RData"Full_Simulation(model_params_true = model_params_true,
true_model = NULL,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 1',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'rule_based',
starting_dose = starting_dose,
SoC = SoC,
use_SoC_data = T)## [1] "Simulation scenario 1_rule_based_well_specified_All_data.RData"toc()## 0.122 sec elapsedcompare_rule_vs_model(sim_title = 'Simulation scenario 1',true_model = NULL,
model_params_true = model_params_true,
prior_model_params = prior_model_params,use_SoC_data = T,
idiosyncratic = F,plot_vstar = T, MTT = MTT, TEL = TEL)
## For the rule-based design, 10% of trials give patient 260 a dose within +/-10% of the true optimal dose
## For the model-based design, 26% of trials give patient 260 a dose within +/-10% of the true optimal doseSimulation 2
true_TED = 8
true_antivenom_mean = 4
true_antivenom_sd = (true_TED - true_antivenom_mean)/qnorm(.95)
true_alpha_tox = logit(1/500) # toxicity at 1 vial
true_MTD = 20 # simulation truth for the MTD
true_beta_tox = solve_beta(alpha_val = true_alpha_tox, v_star = true_MTD, y_star = MTT)
# well-specified simulation truth
model_params_true = list(beta_tox = true_beta_tox,
alpha_tox=true_alpha_tox,
mu_antivenom=true_antivenom_mean,
sd_antivenom=true_antivenom_sd)
optimal_doses[2] = min(true_MTD, true_TED)
tic()
Full_Simulation(model_params_true = model_params_true,
true_model = NULL,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 2',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'model_based',
starting_dose = starting_dose,
SoC = SoC,
use_SoC_data = T)## [1] "Simulation scenario 2_model_based_well_specified_All_data.RData"Full_Simulation(model_params_true = model_params_true,
true_model = NULL,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 2',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'rule_based',
starting_dose = starting_dose,
SoC = SoC,
use_SoC_data = T)## [1] "Simulation scenario 2_rule_based_well_specified_All_data.RData"toc()## 0.063 sec elapsedcompare_rule_vs_model(sim_title = 'Simulation scenario 2',true_model = NULL,
model_params_true = model_params_true,
prior_model_params = prior_model_params,use_SoC_data = T,
idiosyncratic = F,plot_vstar = T, MTT = MTT, TEL = TEL)
## For the rule-based design, 55% of trials give patient 260 a dose within +/-10% of the true optimal dose
## For the model-based design, 83% of trials give patient 260 a dose within +/-10% of the true optimal doseSimulation 3
true_TED = 60
true_antivenom_mean = 30
true_antivenom_sd = (true_TED - true_antivenom_mean)/qnorm(.95)
true_alpha_tox = logit(1/1000) # toxicity at 1 vial
true_MTD = 30 # simulation truth for the MTD
true_beta_tox = solve_beta(alpha_val = true_alpha_tox, v_star = true_MTD, y_star = MTT)
# well-specified simulation truth
model_params_true = list(beta_tox = true_beta_tox,
alpha_tox = true_alpha_tox,
mu_antivenom=true_antivenom_mean,
sd_antivenom=true_antivenom_sd)
optimal_doses[3] = min(true_MTD, true_TED)
tic()
Full_Simulation(model_params_true = model_params_true,
true_model = NULL,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 3',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'model_based',
starting_dose = starting_dose,
SoC = SoC,
use_SoC_data = T)## [1] "Simulation scenario 3_model_based_well_specified_All_data.RData"Full_Simulation(model_params_true = model_params_true,
true_model = NULL,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 3',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'rule_based',
starting_dose = starting_dose,
SoC = SoC,
use_SoC_data = T)## [1] "Simulation scenario 3_rule_based_well_specified_All_data.RData"toc()## 0.081 sec elapsedcompare_rule_vs_model(sim_title = 'Simulation scenario 3',true_model = NULL,
model_params_true = model_params_true,
prior_model_params = prior_model_params,use_SoC_data = T,
idiosyncratic = F,plot_vstar = T, MTT = MTT, TEL = TEL)
## For the rule-based design, 31% of trials give patient 260 a dose within +/-10% of the true optimal dose
## For the model-based design, 27% of trials give patient 260 a dose within +/-10% of the true optimal doseSimulation 4
true_TED = 30
true_antivenom_mean = 15
true_antivenom_sd = (true_TED - true_antivenom_mean)/qnorm(.95)
true_alpha_tox = logit(1/1000) # toxicity at 1 vial
true_MTD = 60 # simulation truth for the MTD
true_beta_tox = solve_beta(alpha_val = true_alpha_tox, v_star = true_MTD, y_star = MTT)
# well-specified simulation truth
model_params_true = list(beta_tox = true_beta_tox,
alpha_tox = true_alpha_tox,
mu_antivenom=true_antivenom_mean,
sd_antivenom=true_antivenom_sd)
optimal_doses[4] = min(true_MTD, true_TED)
tic()
Full_Simulation(model_params_true = model_params_true,
true_model = NULL,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 4',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'model_based',
starting_dose = starting_dose,
SoC = SoC,
use_SoC_data = T)## [1] "Simulation scenario 4_model_based_well_specified_All_data.RData"Full_Simulation(model_params_true = model_params_true,
true_model = NULL,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 4',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'rule_based',
starting_dose = starting_dose,
SoC = SoC,
use_SoC_data = T)## [1] "Simulation scenario 4_rule_based_well_specified_All_data.RData"toc()## 0.065 sec elapsedcompare_rule_vs_model(sim_title = 'Simulation scenario 4',true_model = NULL,
model_params_true = model_params_true,
prior_model_params = prior_model_params,
use_SoC_data = T, idiosyncratic = F,plot_vstar = T,
MTT = MTT, TEL = TEL)
## For the rule-based design, 62% of trials give patient 260 a dose within +/-10% of the true optimal dose
## For the model-based design, 91% of trials give patient 260 a dose within +/-10% of the true optimal doseSimulation 5
Idiosyncratic toxicity
This simulation assumes a dose-independent probability of having a toxic event. We simulate a trial whereby the overall toxicity is 15%
true_model = list(tox = approxfun(x = c(1,10), y = c(.15,.15), rule = 2),
eff = approxfun(x = 0:100,
y = pnorm(q = 0:100, mean = true_antivenom_mean,
sd = true_antivenom_sd),rule=2))
optimal_doses[5] = 0
tic()
# Model based simulation
Full_Simulation(model_params_true = NULL,
true_model = true_model,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 5',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'model_based',
starting_dose = starting_dose,
SoC = SoC,use_SoC_data = T)## [1] "Simulation scenario 5_model_based_mis_specified_All_data.RData"# Rule based simulation
Full_Simulation(model_params_true = NULL,
true_model = true_model,
prior_model_params = prior_model_params,
N_trials = 10*N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 5',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'rule_based',
starting_dose = starting_dose,
SoC = SoC, use_SoC_data = T)## [1] "Simulation scenario 5_rule_based_mis_specified_All_data.RData"toc()## 0.206 sec elapsed# # Comparison in mis-specified idiosyncratic case
compare_rule_vs_model(sim_title = 'Simulation scenario 5',
model_params_true = NULL,
true_model = true_model,
prior_model_params = prior_model_params,
use_SoC_data = T, idiosyncratic = T,plot_vstar = F,
MTT = MTT, TEL = TEL)
## For the rule-based design, 0% of trials give patient 260 a dose within +/-10% of the true optimal dose
## For the model-based design, 0% of trials give patient 260 a dose within +/-10% of the true optimal doseSimulation 6
Same as simulation 2 but with mis-specified efficacy dose response
true_alpha_tox = logit(1/500) # toxicity at 1 vial
true_MTD = 20 # simulation truth for the MTD
true_beta_tox = solve_beta(alpha_val = true_alpha_tox, v_star = true_MTD, y_star = MTT)
true_rate = 1/2.7
# mis-specified simulation truth
true_model = list(eff = approxfun(x =0:100, y = pexp(q = 0:100, rate = true_rate), rule = 2),
tox = approxfun(x = 0:200,
y = inv.logit(true_alpha_tox + true_beta_tox*log2(0:200)),rule=2))
optimal_doses[6] = qexp(p = TEL, rate = true_rate)
tic()
Full_Simulation(model_params_true = NULL,
true_model = true_model,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 6',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'model_based',
starting_dose = starting_dose,
SoC = SoC,
use_SoC_data = T)## [1] "Simulation scenario 6_model_based_mis_specified_All_data.RData"Full_Simulation(model_params_true = NULL,
true_model = true_model,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 6',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'rule_based',
starting_dose = starting_dose,
SoC = SoC,
use_SoC_data = T)## [1] "Simulation scenario 6_rule_based_mis_specified_All_data.RData"toc()## 0.072 sec elapsedcompare_rule_vs_model(sim_title = 'Simulation scenario 6',true_model = true_model,
model_params_true = NULL,
prior_model_params = prior_model_params,use_SoC_data = T,
idiosyncratic = F,plot_vstar = T,MTT = MTT, TEL = TEL)
## For the rule-based design, 37% of trials give patient 260 a dose within +/-10% of the true optimal dose
## For the model-based design, 50% of trials give patient 260 a dose within +/-10% of the true optimal doseSimulation 7
Same as simulation 4 except with exponentially distributed venom volume
true_alpha_tox = logit(1/1000) # toxicity at 1 vial
true_MTD = 60 # simulation truth for the MTD
true_beta_tox = solve_beta(alpha_val = true_alpha_tox, v_star = true_MTD, y_star = MTT)
true_rate = 1/10
# mis-specified simulation truth
true_model = list(eff = approxfun(x =0:200, y = pexp(q = 0:200, rate = true_rate), rule = 2),
tox = approxfun(x = 0:500,
y = inv.logit(true_alpha_tox + true_beta_tox*log2(0:500)),rule=2))
optimal_doses[7] = qexp(p = TEL, rate = true_rate)
tic()
Full_Simulation(model_params_true = NULL,
true_model = true_model,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 7',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'model_based',
starting_dose = starting_dose,
SoC = SoC,
use_SoC_data = T)## [1] "Simulation scenario 7_model_based_mis_specified_All_data.RData"Full_Simulation(model_params_true = NULL,
true_model = true_model,
prior_model_params = prior_model_params,
N_trials = N_trials,
MTT = MTT, TEL = TEL,
N_max = N_max,
N_batch = N_batch,
max_increment = max_increment,
Randomisation_p_SOC = Randomisation_p_SOC,
sim_title = 'Simulation scenario 7',
FORCE_RERUN=FORCE_RERUN,
N_cores = N_cores,
design_type = 'rule_based',
starting_dose = starting_dose,
SoC = SoC,
use_SoC_data = T)## [1] "Simulation scenario 7_rule_based_mis_specified_All_data.RData"toc()## 0.066 sec elapsedcompare_rule_vs_model(sim_title = 'Simulation scenario 7',true_model = true_model,
model_params_true = NULL,
prior_model_params = prior_model_params,
use_SoC_data = T, idiosyncratic = F,plot_vstar = T,
MTT = MTT, TEL = TEL)
## For the rule-based design, 31% of trials give patient 260 a dose within +/-10% of the true optimal dose
## For the model-based design, 25% of trials give patient 260 a dose within +/-10% of the true optimal doseSummary of all trials
specification_type = 'well_specified'
mypallete = brewer.pal(9,'Set1')[c(1,2,7)]
par(las=1,bty='n',family = 'serif', cex.lab=1.5, cex.axis=1.5)
layout(mat = matrix(data = c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,6,6,7,7),nrow = 3,byrow = T))
for(ss in 1:7){
if(ss>4) specification_type = 'mis_specified'
f_name_rule = paste('SimulationOutputs/', 'Simulation scenario ', ss,
'_', 'rule_based', '_', specification_type,
'_', 'All_data', '.RData', sep = '')
f_name_model = paste('SimulationOutputs/','Simulation scenario ', ss,
'_','model_based', '_', specification_type,
'_', 'All_data', '.RData', sep = '')
load(f_name_model)
Summary_model = Summary_trials
load(f_name_rule)
Summary_rule = Summary_trials
rm(Summary_trials)
Summary_model_opt = Summary_model[, grep('assigned_dose', colnames(Summary_model))]
opt_dose = optimal_doses[ss]
Summary_model_opt_qs = 10* (apply(Summary_model_opt,1,quantile,probs = c(0.025,0.975),na.rm=T) - opt_dose)
Summary_rule_opt_qs = 10 * (apply(Summary_rule,1,quantile,probs = c(0.025,0.975),na.rm=T) - opt_dose)
plot(10*(rowMeans(Summary_model_opt,na.rm=T) - opt_dose), type='l', lwd=3,
ylim=range(c(0,Summary_model_opt_qs)),col=mypallete[2],
xlab='Number of patients enrolled',
ylab = 'Difference from optimal (mL)')
polygon(x = c(1:nrow(Summary_model),rev(1:nrow(Summary_model))),
y = c(Summary_model_opt_qs[1,], rev(Summary_model_opt_qs[2,])),
col = adjustcolor(mypallete[2],alpha.f = .2), border = NA)
lines(10*(rowMeans(Summary_rule,na.rm=T) - opt_dose), lwd=3, lty=2,col=mypallete[3])
polygon(x = c(1:nrow(Summary_model),rev(1:nrow(Summary_model))),
y = c(Summary_rule_opt_qs[1,], rev(Summary_rule_opt_qs[2,])),
col = adjustcolor(mypallete[3],alpha.f = .2), border = NA)
mtext(text = ss,side = 3,line = 1,adj = 0)
abline(h=0, col=1)
if(ss==1){
legend('topright',col = mypallete[2:3], legend = c('Model based', 'Rule based'),
cex=1.2,lwd=3,lty=1:2, bty='n')
}
}